Virtual greenhouse model

Niels Holst, Aarhus University, Denmark, 25 May 2023.

R scripts

The documentation comes with R script that demonstrate the equations and logic of the model. Before you run any of the scripts, you need to change the R work folder. This is defined in the setup.R script, which begins like this:

1
# Prelude
2
setwd("~/QDev/UniSim3/input/models/vg")
3
library(ggplot2)
4
library(ggpubr)
5
library(plyr)
6
library(reshape2)
7
...

• Edit the file path in the sewtd(...) statement to point to the folder where you put the R scripts.

• Make sure the four libraries (loaded in lines 3 to 6) are all installed on your computer.

• Run the setup.R script once before running any of the other scripts.

Model structure

The model is updated iteratively in a loop at a time step of $\mathrm{\Delta }t$$\Delta t$ (s), typically in a range of 2 to 5 minutes; however, the energy and water budgets may be resolved at a finer time step as needed (dividing $\mathrm{\Delta }t$$\Delta t$ into $n$$n$ sub-steps). The model runs from the chosen start date until the end date; however, beneath the surface it will always start 24 hours before the start date to allow model state variables to evolve from their initial values and reach realistic values.

Model components are updated in the following order:

1. Calender updates current date and time.

2. Outdoors updates climatic variables based on the weather file. Readings from the weather file are interpolated linearly to estimate values for the current date and time. If the weather file does not contain a Tsky column, sky temperature is estimated by a model. Soil temperature is currently roughly modelled (this needs an update).

3. Greenhouse construction contains the dimensions of the greenhouse and a description of each of the six greenhouse faces (traditionally, glass, see cover layer).

4. Setpoints are defined for the green house actuators, such as growth lights, heat pipes, ventilation, screens and CO₂ diffusers, to keep indoors climate within the desired bounds.

5. Controllers use the setpoints to compute the settings of actuators.

6. Actuators represent growth lights, heat pipes, ventilation, screens, CO₂ diffusers, etc., which control the greenhouse climate. They are governed by the setpoints via the controllers.

7. Indoors holds the greenhouse compartments of layers and volumes and keeps their state variables updated, reflecting the current outdoors driving variables and actuator settings.

8. Output writes the chosen variables to the output file.

Boxscript example with calendar

The simulation begins 24 hours before the declared beginning of the simulation. This allows model state variables to attain sensible initial values. Here is a possible layout of a boxscript. The lonely : characters denote left-out lines:

xxxxxxxxxx
38
1
Simulation sim {
2
  Box global {
3
    &beginDate = 2022/08/01
4
    :
5
  }
6
  Calendar calendar {
7
    .latitude  = 55.4
8
    .longitude = 10.4
9
    .begin = global[beginDate] - 1
10
    .end   = 2022/08/15
11
    .timeStep = 15
12
    .timeUnit = "m"
13
  }
14
  Outdoors outdoors { 
15
  : 
16
  }
17
  Box construction {
18
  :
19
  }
20
  Box setpoints {
21
  :
22
  }
23
  Controllers controllers {
24
  :
25
  }
26
  Box actuators {
27
  :
28
  }
29
  Box indoors {
30
  :
31
  }
32
  OutputR output {
33
    :
34
    OutputSelector {
35
      .beginDateTime = global[beginDate]
36
    }
37
  }
38
}

A global box is declared to hold global settings, i.e., parameters that serve as inputs to one or more of the boxes further down in the boxscript. A beginDate port is created (the ampersand & means that a new port is created) to start the simulation on 1 August 2022. The calendar actually starts one day before, as set by its begin port (line 9); the initial period . denotes an existing input port. At the very end of the boxscript, an OutputSelector box asks for the output to be written starting on 1 August (line 35). In effect, the initial extra day's output will be discarded. The simulation will run until the final date end set by the calendar box; the time step has been set to 15 minutes (lines 10-11).

Greenhouse virtualization

Greenhouse compartments

Greenhouse dimensions

The greenhouse consists of $n$$n$ spans each of width $w$$w$ (m), length $\ell$$\ell$ (m) and wall height $h$$h$ (m). Each span has a roof on top angled at $u$$u$ (°).

From these parameters we can derive several measures of greenhouse geometry. The area of a single face, either side, end or roof, is

from which follows the total area of side, end and roof faces, and the total cover area (all in m²),

The ground area (m²) enclosed by the greenhouse is

and the volume contained (m³) is

Further useful measures are the cover/area ratio $C_{ai}$$C_{ai}$ (m² cover/m² ground),

and the volumetric area index $V_{ai}$$V_{ai}$ (m) (also known as the average height),

Furthermore, the density of air ${\rho }_{air}=1.19$$\rho_{air}=1.19$ kg/m³ and it heat capacity $C_{air}=1020$$C_{air}=1\,020$ J/kg/K leads to the heat capacity of indoors air $C_{gh}$$C_{gh}$ (J/K/m² ground),

The equations above have been coded in the greenhouse-geometry.R script::

​x
1
n = 25
2
w = 4
3
L = 100
4
h = 3.5
5
u = 26
6

7
A_side_            350.000
8
A_end_              15.753
9
A_roof_            218.373
10

11
A_side             700.000
12
A_end              787.674
13
A_roof           10918.650
14
A_cover          12406.325
15

16
A_gh             10000.000
17
V_gh             39383.711
18
cover_per_ground     1.241
19
avg_height           3.938
20
C_gh              4780.395

Boxscript example

In boxscripts, greenhouse dimensions are set in the geometry box inside construction:

xxxxxxxxxx
10
1
Box construction{
2
  vg::Geometry geometry{
3
    .numSpans = 25
4
    .spanWidth = 4
5
    .length = 100
6
    .height = 3.5
7
    .roofPitch = 26
8
  }
9
  :
10
}

Volumes

The greenhouse model holds three volumes: indoors, outdoors and soil. The indoors holds state variables, while the outdoors and soil hold driving variables. All three hold temperature $T$$T$ (°C), while indoors and outdoors hold relative humidity $H^{\mathrm{\%}}$$H^\%$ (%) and CO₂ concentration $[C{O}_2{]}^{ppm}$$[CO_2]^{ppm}$ as well. In some relations other units are more convenient, e.g. absolute humidity $H^{abs}$$H^{abs}$ (kg/m³), water vapour pressure $H^{pres}$$H^{pres}$ (Pa), or absolute CO₂ concentration $[C{O}_2{]}^{abs}$$[CO_2]^{abs}$ (g/m³). The variable are subscripted to denote the volume, e.g., for temperature $T_{in}$$T_{in}$, $T_{out}$$T_{out}$ and $T_{soil}$$T_{soil}$.

Volumes may exchange heat and humidity with neighbouring volumes and layers.

Layers

Besides layers the greenhouse and its environment is conceptualised as a stack of $i=1..n$$i=1..n$ horizontal layers. This is a typical stack:

1. sky

2. cover (e.g., glass)

3. screen

4. growth light

5. plant canopy

6. heat pipe

7. floor

Depending on the scenario, layers can be removed or added (e.g., additional screen layers or a table layer ). All layers, whether they are horisontal or inclined, uniform or heterogeneous, thick or thin, are simplified to describe an average 1m² ground area of the greenhouse. For screens this area will vary according to what degree they are drawn. Where the unit m² is used, it always refers to m² ground. Sometimes it is spelled out "m² ground" to avoid confusion.

The radiative properties of each layer is parameterised for three intervals of the spectrum:

• PAR (μmol/m²/s)

• short-wave radiation (W/m²)

• long-wave radiation (W/m²)

PAR is a subset of the short-wave radiation. The PAR budget is resolved only to compute how much is absorbed by the plant canopy. Its quantum units reflect this purpose. In contrast, the short- and long-wave radiation both contribute to the radiative exchange of heat among layers. Remember that W = J/s.

Each layer $i$$i$ is characterised by its radiative properties: absorptivity $({\alpha }_i^r$$(\alpha^r_i$), reflectivity $({\rho }_i^r)$$(\rho^r_i)$ and transmissivity $({\tau }_i^r)$$(\tau^r_i)$ for radiation coming from above and similarly for radiation coming from below $(\alpha {{}_i^r}^{\prime },\rho {{}_i^r}^{\prime },\tau {{}_i^r}^{\prime })$$(\alpha{^r_i}',\rho{^r_i}',\tau{^r_i}')$, all for each type of radiation $r$$r$ (PAR, sw, lw). The radiative parameters always add up to unity, ${\alpha }_i^r+{\rho }_i^r+{\tau }_i^r=1$$\alpha^r_i+\rho^r_i+\tau^r_i=1$ and $\alpha {{}_i^r}^{\prime }+\rho {{}_i^r}^{\prime }+\tau {{}_i^r}^{\prime }=1$$\alpha{^r_i}'+\rho{^r_i}'+\tau{^r_i}'=1$. We apply the common engineering assumption that emissivity equals the absorptivity for long-wave radiation, i.e. for all layers ${ϵ}_i={\alpha }_i^{lw}$$\epsilon_i=\alpha_i^{lw}$ and ${ϵ}_i^{\prime }={{\alpha }_i^{lw}}^{\prime }$$\epsilon_i'={\alpha_i^{lw}}'$.

The superscript $r$$r$ was left out of the figure above for simplification but it is important, as it allows different parameter settings for different radiation intervals. This is essential, e.g. for the distinct difference in transmissivity of glass for short- vs. long-wave radiation. Moreover, the separate parameter sets for down- and upgoing radiation allow for asymmetry, e.g., of screens.

Layers are further described by their area-specific heat capacity ($C_i$$C_i$; J/K/m² ground) and temperature ($T_i$$T_i$; °C).

Boxscript example with cover and screens

The construction box holds two child boxes, geometry and shelter. The shelter box holds definitions of different cover and screen materials inside the covers and screens box, respectively. In the boxscript below, one kind of cover is defined (named glass) and four kinds of screens (named energy, shade1, shade2 and blackout). You can name the materials as you want.

The faces box (lines 67-98) holds the named materials of cover and screens for each of the six shelter faces. There may be more than one layer of screens installed. These screen layers are written in sequence starting nearest the glass. For example, with up to three layers installed on a face, you might have faces with screens specified as energy+shade1+blackout, none+shade2+blackout, energy+none+blackout or none+none+blackout, or any other combination of screen materials defined by the names of the shelter/screens/* boxes. All six faces must have the same number of screens defined. If a face has no screen in a certain layer, it is written as none.

xxxxxxxxxx
106
1
Box construction{
2
  vg::Geometry geometry{
3
    .numSpans = 25
4
    .spanWidth = 4
5
    .length = 100
6
    .height = 3.5
7
    .roofPitch = 26
8
  }
9
  Box shelter {
10
    UWind Utop {
11
    }
12
    Box covers {
13
      Cover glass {
14
        .swReflectivityTop      = 0.1
15
        .swTransmissivityTop    = 0.8
16
        .swReflectivityBottom   = 0.1
17
        .swTransmissivityBottom = 0.8
18
        .lwReflectivityTop      = 0.15
19
        .lwTransmissivityTop    = 0.20
20
        .lwReflectivityBottom   = 0.15
21
        .lwTransmissivityBottom = 0.20
22
        .Utop                   = ../../Utop[value]
23
        .Ubottom                = 1.2
24
        .heatCapacity           = 8400
25
        .swReflectivityChalk    = actuators/chalk[swReflectivity]
26
        .lwReflectivityChalk    = actuators/chalk[lwReflectivity]
27
      }
28
    }
29
    Box screens {
30
      Screen energy {
31
        .swReflectivityTop      = 0.93
32
        .swTransmissivityTop    = 0.0
33
        :
34
        .Utop                   = 1.2
35
        .Ubottom                = 1.2
36
        .heatCapacity           = 80
37
        .state = actuators/screens/energy[state]
38
      }
39
      Screen shade1 {
40
        .swReflectivityTop      = 0.70
41
        .swTransmissivityTop    = 0.20
42
        :
43
        .Utop                   = 1.2
44
        .Ubottom                = 1.2
45
        .heatCapacity           = 80.
46
        .state = actuators/screens/shade[state]
47
      }
48
      Screen shade2 {
49
        .swReflectivityTop      = 0.80
50
        .swTransmissivityTop    = 0.10
51
        :
52
        .Utop                   = 1.1
53
        .Ubottom                = 1.1
54
        .heatCapacity           = 120.
55
        .state = actuators/screens/shade[state]
56
      }
57
      Screen blackout {
58
        .swReflectivityTop      = 0.05
59
        .swTransmissivityTop    = 0.0
60
        :
61
        .Utop                   = 1.2
62
        .Ubottom                = 1.2
63
        .heatCapacity           = 200.
64
        .state = actuators/screens/blackout[state]
65
      }
66
    }
67
    Box faces {
68
      Box roof1 {
69
        &cover = "glass"
70
        &screens = "energy+shade1+blackout"
71
        &area = ../../geometry[roofArea]/2
72
        &weight = 1
73
      }
74
      Box roof2 {
75
        &cover = "glass"
76
        &screens = "energy+shade1+blackout"
77
        &area = ../../geometry[roofArea]/2
78
        &weight = 1
79
      }
80
      Box side1 {
81
        &cover = "glass"
82
        &screens = "none+shade2+blackout"
83
        &area = ../../geometry[sideArea]/2
84
        &weight = 1
85
      }
86
      Box side2 {
87
        &cover = "glass"
88
        &screens = "none+none+blackout"
89
        &area = ../../geometry[sideArea]/2
90
        &weight = 0.5
91
      }
92
      Box end1 {
93
        &cover = "glass"
94
        &screens = "none+none+blackout"
95
        &area = ../../geometry[endArea]/2
96
        &weight = 0.2
97
      }
98
      Box end2 {
99
        &cover = "glass"
100
        &screens = "energy+none+blackout"
101
        &area = ../../geometry[endArea]/2
102
        &weight = 0.2
103
      }
104
    }
105
  }
106
}

Any layer has 11 input ports defining its basic radiative properties (lines 14-21), $U$$U$-values (lines 22-23) and heat capacity (line 24). Cover boxes holds two additional inputs (lines 25-26) specifying the current status of whitening ('chalk'), which is provided by the chalk actuator. In a similar vein, Screen boxes hold one additional input state, which refers to a screen actuator. From the code above, you can infer that in this model there are three screen actuators:

• actuators/screens/energywhich determines the position of the energy screen

• actuators/screens/shade which determines the position of the shade1 and shade2 screens

• actuators/screens/blackout which determines the position of the blackout screen

In addition to cover and screens, each face has a weight parameter that is used when calculating the average radiative properties of the greenhouse cover and each screen layer. The settings above suggest that side2 is north-facing.

Processes of heat exchange

Radiative exchange

Radiation might be emitted by a layer, both down- and upwards, in terms of PAR ($E_i^{par},{E_i^{par}}^{\prime }$$E_i^{par},\,{E_i^{par}}'$; μmol/m²/s) and short- and long-wave radiation ($E_i^{sw},{E_i^{sw}}^{\prime },E_i^{lw},{E_i^{lw}}^{\prime }$$E_i^{sw},\,{E_i^{sw}}',\, E_i^{lw},\,{E_i^{lw}}'$; W/m²). For some layers (sky, cover, screen, floor) , the emission of long-wave radiation is calculated from their temperature, according to the Stefan-Boltzmann Law,

where $T_0=273.15$$T_0=273.15$ °C converts to absolute temperature. For other layers (growth light, plant canopy, heat pipe), the emission is calculated by specific sub-models..

The emitted radiation is intercepted by neighbouring layers, its ultimate fate determined by the characteristics (absorptivity, reflectivity, transmissivity) of the whole stack of layers, as worked out in the energy budget. The net radiative absorption of a layer (i.e. its absorption minus emission) is denoted by ($A_i^{par},{A_i^{par}}^{\prime }$$A_i^{par},\,{A_i^{par}}'$; μmol/m²/s) and ($A_i^{sw},{A_i^{sw}}^{\prime },A_i^{lw},{A_i^{lw}}^{\prime }$$A_i^{sw},\,{A_i^{sw}}',\, A_i^{lw},\,{A_i^{lw}}'$; W/m²).

Convection and conduction

The exchange of heat between layer $i$$i$'s upper and lower surface and the confronting volume is described by $U_i$$U_i$ on the upper side and $U_i^{\prime }$$U_{i}'$ on the lower side (both in W/K/m²). The mechanism of heat exchange is convection for air volumes and conduction for the soil volume but despite that they are all superscripted as $A_i^{conv}$$A_i^{conv}$.

The convective/conductive heat absorbed by layer $i$$i$ (W/m²) from the volume above and below is, respectively,

where the temperature of the volume above and below layer $i$$i$ is denoted by $T_{Vi}$$T_{Vi}$ and $T_{Vi}^{\prime }$$T_{Vi}'$, respectively.

Advection

Ventilation rate integration

Advective heat transfer occurs when outdoors air replaces indoors air. This is due to intentional ventilation as well as inadvertent leakage. We formulate the ventilation rate as, "how many times greenhouse air is replaced per hour" $v$$v$ (h${}^{-1}$$^{-1}$).

During a time step $\mathrm{\Delta }t$$\Delta t$ the indoors temperature $T_{in}$$T_{in}$ and absolute humidity $H_{in}^{abs}$$H_{in}^{abs}$ will approach the outdoors $T_{out}$$T_{out}$ and $H_{out}^{abs}$$H_{out}^{abs}$. Considering either variable and assuming a perfect mixture of air inside the greenhouse, the rate of change in the indoors variable $y_{in}$$y_{in}$ is

Integration gives us

The relation is illustrated by the advection.R script with $T_{out}$$T_{out}$ at 20°C and $T_{in}(0)$$T_{in}(0)$ either 15°C or 25°C. The ventilation rate was set at $v=3$$v=3$ ${\text{h}}^{-1}$$\text{h}^{-1}$.

To calculate the change in indoors temperature $\mathrm{\Delta }{T}_{in}$$\Delta T_{in}$ during a time step $\mathrm{\Delta }t$$\Delta t$ (s) eq. $\text{14}$$\ref{eq_advection}$ can be rewritten as

With a time step of, for example $\mathrm{\Delta }t=2$$\Delta t=2$ minutes, the decrease (blue curve) or increase (red curve) in indoors temperature would be $\mathrm{\Delta }{T}_{in}=0.48$$\Delta T_{in}=0.48$ K (by eq. $\text{15}$$\ref{eq_advection_temp}$). Assuming a volumetric area index of $V_{ai}=3.94$$V_{ai}=3.94$ m³/m² (eq. $\text{10}$$\ref{average_height}$) and with air heat capacity $C_{air}=1214$$C_{air}=1\,214$ J/m³/K, we get the following energy flux (positive in the case of the red curve),

The change in indoors absolute humidity $\mathrm{\Delta }{H}_{in}^{abs}$$\Delta H_{in}^{abs}$ can be arrived at in a similar fashion.

Eq. $\text{15}$$\ref{eq_advection_temp}$ is an example of a local integral; it is solved assuming everything else remains constant. It is applied here because high values of $v$$v$ could make a straightforward Euler integration imprecise. In this case, the difference is not so great. Euler integration would have yielded

Maximum ventilation rate

The maximum ventilation rate $v_{max}$$v_{max}$ (h${}^{-1}$$^{-1}$) is obtained when the vents are fully open. What this ventilation rate amounts to depends on the wind speed $w$$w$ (m/s) and the difference between the indoors $T_{in}$$T_{in}$ and outdoors $T_{out}$$T_{out}$ temperature (°C):

with discharge coefficients for wind $c_w$$c_w$ (s/m/h) and temperature $c_T$$c_T$ (/K/h). The total discharge is found by vector addition, considering the two forces operating independently (orthogonally). It is applied to the relative vent area; $A_{vent}$$A_{vent}$ (m²) is the total area of vents and $A_{gh}$$A_{gh}$ (m²) is the greenhouse floor area.

Latent heat

Latent heat is related to the phase change between water in the gas and the liquid phase. The phase change costs (by evaporation) or releases (by condensation) energy according to the latent heat of water vapourisation $\lambda =2454$$\lambda=2\,454$ kJ/kg.

In the model, transpiration cools the canopy and condensation heats the cover. The transpired water adds to air humidity, which is then considered to have increased its content of latent heat. Water vapour bears latent heat in the sense that if it condenses, the heat used for the initial phase change from liquid to gas, will be used in the phase change back from gas to liquid.

The model does not include the cooling effect of water evaporating from wet surfaces, nor the wetting and drying of screens.

Layer characterization

Sky layer

The radiation emitted downwards from the sky includes $E_{sky}^{par}$$E_{sky}^{par}$, $E_{sky}^{sw}$$E_{sky}^{sw}$ and $E_{sky}^{lw}$$E_{sky}^{lw}$. PAR and short wave radiation originates from the sun, while long wave radiation is emitted from the sky as such.

While $E_{sky}^{sw}$$E_{sky}^{sw}$ is measured routinely, this is not the case for $E_{sky}^{par}$$E_{sky}^{par}$. It can be approximated by $E_{sky}^{par}=0.5E_{sky}^{sw}$$E_{sky}^{par}=0.5E_{sky}^{sw}$.

Sky emissivity ${ϵ}_{sky}$$\epsilon_{sky}$ (and thereby sky absorptivity ${\alpha }_{sky}$$\alpha_{sky}$) is estimated from the dew point temperature $T_{dew}$$T_{dew}$ (°C) by the empirical relation (ref.)

Sky temperature $T_{sky}$$T_{sky}$ (°C) is then derived from outdoors temperature $T_{out}$$T_{out}$ (°C) according to Karn et al. 2019

where $T_0=273.15$$T_0=273.15$ °C.

Sky reflectivity ${\rho }_{sky}$$\rho_{sky}$ is zero and consequently transmissivity is ${\tau }_{sky}=1-{ϵ}_{sky}$$\tau_{sky}=1-\epsilon_{sky}$.

Cover layer

The cover layer is commonly made out of glass with the following typical characteristics:

In general, the six faces of the greenhouse may carry different covers with different radiative properties. For the cover layer as a whole, each radiative parameter is summarised as the average over all faces, weighed by face area.

The density of window glass is typically $2.5$$2.5$ kg/m²/mm, which with a thickness of $4$$4$ mm gives a density of $10$$10$ kg/m² cover. The heat capacity of glass is $840$$840$ J/kg/K (Engineering Toolbox, 2023). If we furthermore assume a cover/ground ratio of $1.24$$1.24$ m² cover/m² ground (eq. $\text{9}$$\ref{eq_cai}$), we get a heat capacity of the cover of

The cover emits long-wave radiation symmetrically down and up according to the Stefan-Boltzmann Law (eq. $\text{12}$$\ref{stefan_boltzmann}$).

The $U$$U$-values over both the inner and outer surface can be expected to increase with the wind; immobile air is a very good insulator. For the outer surface of the cover, we use the empirical relation,

where $w$$w$ is the outdoors wind speed (m/s). The $U$$U$-value of the inner surface is set to $1.2$$1.2$ W/K/m².

Cover averaging

Since the six cover faces may differ in their characteristics, we may need to average their parameter values. Some examples applying the equations below of are shown in the cover-averaging.R script.

Cover averaging of radiative parameters

Let's assume that the two end faces differ from the other four faces; we set $\alpha =0.15$$\alpha=0.15$, $\rho =0.35$$\rho=0.35$ and $\tau =0.60$$\tau=0.60$ at the ends walls, while elsewhere, $\alpha =0.10$$\alpha=0.10$, $\rho =0.10$$\rho=0.10$ and $\tau =0.80$$\tau=0.80$. In the full model, these values have been corrected for the effects of chalking before the averaging. Furthermore, let's weigh the end walls by $w_{end}=0.25$$w_{end}=0.25$ and one side by $w_{side1}=0.20$$w_{side1}=0.20$ while the rest have full weight $w_{other}=1.0$$w_{other}=1.0$.

For the radiative parameters ($\alpha$$\alpha$, $\rho$$\rho$ and $\tau$$\tau$), we get their weighted averages (${\alpha }_{cover}$$\alpha_{cover}$, ${\rho }_{cover}$$\rho_{cover}$ and ${\tau }_{cover}$$\tau_{cover}$) for the cover as whole by

where subscript $i$$i$ stands for the six cover faces, and $A_i$$A_i$ (m²) is the area of each face. Since ${\alpha }^0+{\rho }^0+{\tau }^0=1$$\alpha^0+\rho^0+\tau^0=1$ is not guaranteed, these values are subsequently rescaled to a sum of $1$$1$.

With the numbers above and the illustrative greenhouse dimensions we have (with parameters pair-wise in the order: roofs, sides, ends)

which gives us the averaged values,

The roof faces are dominating by their size; hence the average parameter values are close to the roof values. If we change all weights to one $w_i=1$$w_i=1$, we get a larger impact of the ends: ${\alpha }_{cover}=0.102$$\alpha_{cover}=0.102$, ${\rho }_{cover}=0.115$$\rho_{cover}=0.115$ and ${\alpha }_{cover}=0.783$$\alpha_{cover}=0.783$.

Cover averaging of heat capacity

Heat capacity is averaged to apply to the greenhouse ground area (eq. $\text{18}$$\ref{eq_c_glass}$). If, for example, the largest area is covered by glass with a heat capacity of $C_{glass}=8400$$C_{glass}=8\,400$ J/K/m² cover, while the ends have a more sturdy material with $C_{sturdy}=18000$$C_{sturdy}=18\,000$ J/K/m² cover, we can calculate the total heat capacity of the cover per greenhouse area $C_{cover}$$C_{cover}$ (J/K/m² ground) as

where $A_{gh}$$A_{gh}$ (m²) is the ground area of the greenhouse. We apply no weights to the different greenhouse faces, other than their areas, when averaging the heat capacity.

Using the the illustrative greenhouse dimensions we have

which gives us

As expected, this heat capacity is larger than the one calculated for and all-glass greenhouse $C_{cover}=10421$$C_{cover}=10\,421$ J/K/m² ground in eq. $\text{18}$$\ref{eq_c_glass}$.

Cover averaging of U-values

The $U_i$$U_i$ values (W/K/m² cover) are conductances (one for each face, $i=1..6$$i=1..6$). They can also be looked upon as resistances by taking their inverse, $R_i=1/U_i$$R_i=1/U_i$. If we assume that the six resistances work like electrical resistors in parallel, we get the total resistance of the greenhouse cover $R_{cover}$$R_{cover}$ as

It follows that conductances are simply additive. Thus the conductance of the greenhouse cover as a whole is arrived at by the sum of the faces. Relating this to the ground area of the greenhouse $A_{gh}$$A_{gh}$ (m² ground), we get

As for heat capacity (eq. $\text{21}$$\ref{eq_c_cover}$), we apply no weights to the different greenhouse faces, other than their areas.

For example, If the cover material has $U_{glass}=1.2$$U_{glass}=1.2$ W/K/m², except for the end walls which have $U_{poly}=0.3$$U_{poly}=0.3$ W/K/m², we get with the illustrative greenhouse dimensions,

from which follows

The reason why the average $U_{cover}$$U_{cover}$ is larger than any $U_i$$U_i$ is that the ground area is smaller than the cover area.

Screen layer

Screens are often constructed with asymmetric radiative properties. Here, for example, are the radiative parameters of a polyester screen with an aluminium surface on the upper side (from Bailey 1981):

The values in the table are for long-wave radiation only. They are most likely different for short-wave radiation but documentation for the physical characteristics of screens are hard to come by; manufacturers seem reluctant to provide the information.

Screens may be more or less drawn . If we denote how much a curtain is drawn by its position $p\in [0;1]$$p\in[0;1]$ then the effective radiation parameters are calculates as $p\alpha$$p\alpha$, $p\rho$$p\rho$ and $1-p\alpha -p\rho$$1-p\alpha-p\rho$, and likewise for ${\alpha }^{\prime }$$\alpha'$, ${\rho }^{\prime }$$\rho'$ and ${\tau }^{\prime }$$\tau'$. For the screen layer of the greenhouse as a whole, each radiative parameter is summarised as the average over all the six faces, weighed by face area, as for the cover layer. If screens are kept half-drawn for extended periods, or if faces differ much in screen material, this averaging could turn out highly unrealistic. If the greenhouse is equipped with more than one layer of screens, the screens are stacked as separate, independent layers.

Screens may be constructed of composite materials which complicates the calculation of their heat capacity. For the screen exemplified here, we found quite close values of $C_{polyester}=1.1$$C_{polyester}=1.1$ J/K/g and $C_{Al}=0.89$$C_{Al}=0.89$ J/K/g. Assuming that polyester is the dominant material and a screen density of $80$$80$ g/m² screen, we get a heat capacity of the screen of

This calculation implies that all six greenhouses faces are equally covered with the same screen. If that is not the case then $S_{screen}$$S_{screen}$ must be found by calculating the heat capacity of each face separately and adding them up to give $S_{screen}$$S_{screen}$. Note that the calculated heat capacity of the screen is ~100 times less that of a glass cover (eq. $\text{18}$$\ref{eq_c_glass}$).

The screen emits long-wave radiation according to the Stefan-Boltzmann Law (eq. $\text{12}$$\ref{stefan_boltzmann}$). With the example above, the radiation would be asymmetric because downwards and upwards emissivity differ; this screen was designed to keep the heat indoors.

Screen averaging

The averaging of screens in a layer follows the same principles as cover averaging, except that screens also have a position $p\in [0;1]$$p\in[0;1]$. In a given layer, there may be screens on some layers but not on others; these are given full transparency, $\tau ={\tau }^{\prime }=1$$\tau=\tau'=1$.

Thus for the radiative parameter $q$$q$ (i.e. $\alpha$$\alpha$, $\rho$$\rho$, ${\alpha }^{\prime }$$\alpha'$ or ${\rho }^{\prime }$$\rho'$), we get the weighted average $q_{screen}$$q_{screen}$ for the screen layer as a whole (cf. eq. $\text{20}$$\ref{eq_q_cover}$),

with transmissivity as in eq. $\text{???}$$\ref{eq_tau_cover}$,

The heat capacity and U-value of screens is assumed to be in effect only when they are drawn. Hence we get (cf. eqs. $\text{21}$$\ref{eq_c_cover}$ and $\text{22}$$\ref{eq_u_cover}$) the total heat capacity of the screen layer $C_{screen}$$C_{screen}$ (J/K/m² ground),

and the total U-value of the screen layer $U_{screen}$$U_{screen}$ (W/K/m² ground),

Some examples are shown in the screen-averaging.R script.

Screen layering

Screen layers function just like any other layer (see Layers). Thus there is need for special mechanisms to account for screen layering.

Growth light layer

Lamps are installed to provide growth light. Even so, their production of heat is a significant side effect. Light is traditionally emitted downwards from a position above the canopy. If the installed lamp power (including any ballast or driver needed by the lamps) is $P_{lamp}$$P_{lamp}$ (W/m²) and the PAR efficiency is $e_{lamp}$$e_{lamp}$ (μmol PAR/J) then the emission of PAR downwards $E_{lamp}^{par}$$E_{lamp}^{par}$ and upwards ${E_{lamp}^{par}}^{\prime }$${E_{lamp}^{par}}'$ (μmol PAR/m²) is

Lamps may loose efficiency with age, which is reflected in a reduced $e_{lamp}$$e_{lamp}$.

Energy is dissipated from the lamps by three routes: short-waved and long-waved radiation and convection. The proportions are denoted $p_{lamp}^{sw}$$p_{lamp}^{sw}$, $p_{lamp}^{lw}$$p_{lamp}^{lw}$ and $p_{lamp}^{conv}$$p_{lamp}^{conv}$, where $p_{lamp}^{sw}+p_{lamp}^{lw}+p_{lamp}^{conv}=1$$p_{lamp}^{sw}+p_{lamp}^{lw}+p_{lamp}^{conv}=1$. While the short-wave radiation is assumed to be heading only downwards, the long-wave radiation and convective heat are both assumed to be emitted equally downward and upwards. Thus we have

The growth light layer has a 100% transmissivity ($\tau ={\tau }^{\prime }=1$$\tau=\tau'=1$) and is attributed neither a temperature nor a heat capacity.

Plant canopy layer

Stanghellini (1987 ,eq. 2.1) sets up the following energy balance for the plant canopy

The equation states (in her notation) that the energy of net radiation absorbed by the canopy $R_n$$R_n$ is spent on

• convective heat loss to the greenhouse air $H$$H$

• heat for the evaporation of water needed for transpiration $LE$$LE$

• photosynthate production $M$$M$

• heat storage in the canopy $J$$J$

Like Stanghellini we will ignore $M$$M$ and $J$$J$ in the energy budget (they both make quite small contributions to the total budget) and follow her logic which leads to, surprisingly, that $H$$H$ does not need to be taken explicitly into account either.

Radiative properties

Goudriaan (1977) provides equations (his eqs. 2.21 and 2.26) to calculate the reflectivity of a plant canopy with leaf area index $L_{ai}$$L_{ai}$ and extinction coefficient $k$$k$. The reflection coefficient ${\rho }_h$$\rho_h$ depends on the scattering coefficient $\sigma$$\sigma$,

$\sigma$$\sigma$ is a species-dependent parameter commonly set to $\sigma =0.2$$\sigma=0.2$ (Kropff and Laar 1993, p.38), which gives ${\rho }_h=0.0557$$\rho_h=0.0557$.

The original eq. 2.26 includes the reflectivity of the surface underneath the canopy (the floor in our case). However, we will resolve the distribution of radiation among the greenhouse layers otherwise. Consequently, we can set the reflectivity of the underlying surface in the original equation to zero. Thus we get the reflectivity of the canopy, which is symmetric for up- and down-going radiation, based on his eq.2.26,

where superscript $r$$r$ denotes $par$$par$, $sw$$sw$ or $lw$$lw$ radiation. Commonly found values for $k^{par}$$k^{par}$ and $k^{sw}$$k^{sw}$ are in the range $0.7$$0.7$ to $0.8$$0.8$, while leaves are impenetrable to long-wave radiation resulting in $k^{lw}=1$$k^{lw}=1$.

The absorptivity of the canopy is simply

and the transmissivity

The plant-canopy-layer.R script shows how the parameters change with increasing leaf area index for short-wave ($k=0.7$$k=0.7$) and long-wave radiation ($k=1$$k=1$):

Canopy temperature

Note: The figures in this sub-section were generated by the 5-energy-budget-plant.R script.

The model of canopy temperature is taken from Stanghellini (1987 ,eq. 3.5):

with

• plant and indoors climate variables

  • $T_{plant}$$T_{plant}$ : Canopy temperature (°C)

  • $T_{in}$$T_{in}$ : Indoors air temperature (°C)

  • $L_{ai}$$L_{ai}$ : Leaf area index (m² leaf/m² ground)

  • $r_i=r_i^{H_{2}O}$$r_i=r_i^{H_2O}$ : Internal (stomatal) leaf resistance against water vapour (s/m)

  • $r_e=r_e^{H_{2}O}$$r_e=r_e^{H_2O}$ : External (boundary layer) leaf resistance against water vapour (s/m)

  • $A_{plant}^{rad}$$A_{plant}^{rad}$: Net radiation absorbed by the canopy (W/m²)

  • $H_{in}^{pres}$$H_{in}^{pres}$ : Indoors air vapour pressure (Pa)

  • $H_{sat}^{pres}$$H_{sat}^{pres}$ : Saturated air vapour pressure (Pa), a function of $T_{in}$$T_{in}$

• physical variables

  • $\delta$$\delta$ : Slope of saturated vapour pressure curve (Pa/K)

• physical constants and approximations

  • $\lambda$$\lambda$: Latent heat of vapourisation of water ($2454$$2\,454$ kJ/kg )

  • $\gamma$$\gamma$: The psychrometric constant, which tells how the partial pressure of water vapour changes with changes in air temperature ($67.71$$67.71$ Pa/K at standard air pressure)

  • ${\rho }_a$$\rho_a$: Air density, which depends on air pressure, air humidity and temperature (Shelquist, 2023); here assumed constant ($1.19$$1.19$ kg/m³)

  • $C_a$$C_a$: Air heat capacity, which depends on air humidity and temperature (Powder Process, 2023); here assumed constant ($1020$$1\,020$ J/kg/K)

The saturated vapour pressure $H_{sat}(T)$$H_{sat}(T)$ as a function of temperature is described by the empirical equation (FAO, 2023, eq. 11),

It increases with temperature, as seen in the figure:

The slope $\delta (T)$$\delta(T)$ as a function of temperature is described by the empirical equation (FAO, 2023, eq. 13),

The slope is increasing with temperature too, as seen in the figure:

The internal leaf (stomatal) resistance against water vapour $r_i^{H_{2}O}$$r_i^{H_2O}$ is computed by the Ball-Barry (ref.) function,

with species-specific coefficients that we here set to $g_0=0.1$$g_0=0.1$ m/s and $g_1=1.64$$g_1=1.64$ m³/mol. Other variables are

• $H_{in}^{rel}$$H_{in}^{rel}$: Relative humidity of the indoors air (%)

• $c_{in}^{C{O}_2}$$c_{in}^{CO_2}$: Indoors CO₂ concentration (ppm)

• $P_n:$$P_n:$ Net photosynthetic rate (μmol CO₂ / leaf m² / s)

Here is an example, how $r_i^{H_{2}O}$$r_i^{H_2O}$ behaves with increasing humidity at various CO₂ levels and with $P_n=2.0$$P_n=2.0$ μmol CO₂ / m² leaf / s. There is not much variation from the maximum set by $g_0^{-1}$$g_0^{-1}$:

The external leaf (boundary layer) resistance to water vapour $r_e^{H_{2}O}$$r_e^{H_2O}$ can be computed according to Stanghellini (1987, eq. 2.53):

where $\ell$$\ell$ characterizes leaf dimensions and $u$$u$ (m/s) is air velocity at the leaf surface. The constants $1174$$1\,174$ and $207$$207$ were estimated empirically in a tomato crop by Stanghellini.

This figure shows the behaviour of $r_e^{H_{2}O}$$r_e^{H_2O}$ with increasing wind speed, at different leaf temperatures and leaf dimensions with $T_{in}=24$$T_{in}=24$°C:

Both $\ell$$\ell$ and $u$$u$ are difficult to estimate. $\ell$$\ell$ should be smaller for smaller leaves. $u$$u$ should increase with the greenhouse ventilation rate (or if rotors are running) and decrease with $L_{ai}$$L_{ai}$, as a denser canopy will reduce the average wind speed at leaf surfaces.

Note: Due to all the uncertainties, we use a fixed value in the model defaulting to

Everything above taken together, we get the following linear response of plant temperature $T_{plant}$$T_{plant}$ (eq. $\text{27}$$\ref{eq_t_plant}$) to the radiation absorbed by canopy $A_{plant}^{rad}$$A_{plant}^{rad}$ with $L_{ai}=1.9$$L_{ai}=1.9$. The indoors temperature was set at $T_{in}=25$$T_{in}=25$°C and is shown as a black line in the figure:

Clearly, a high air humidity works against plant transpiration and leads to a higher leaf temperature, even above the ambient temperature of the greenhouse: Relatively cool leaves are a sign of a sound microclimate.

Transpiration

Note: The figures in this sub-section were generated by the 5-energy-budget-plant.R script.

The model of canopy transpiration rate $H_{plant}^{trans}$$H_{plant}^{trans}$ (kg/m² ground/s) is taken from Stanghellini (1987 ,eq. 3.4):

with symbols defined as for eq. $\text{27}$$\ref{eq_t_plant}$. With the same parameter settings as for the previous figure, we get this transpiration rate of the canopy (shown in units of g/m² ground/min):

Obviously, transpiration is larger in a drier climate and it increases with the absorbed radiation. Note that leaf temperature (previous figure) too increased with absorbed radiation. This means that the increase in transpiration with increased absorbed radiation was not sufficient to keep leaf temperature constant; that would have acquired an even higher transpiration.

The phase change from water (in the plant) to water vapour (in the greenhouse air) causes an increase in the latent heat carried by the air. This heat, which is equal to $\lambda H_{plant}^{trans}$$\lambda H_{plant}^{trans}$, is shown in the figure below,

The 1:1 dashed line has been added to highlight the balance between the radiation absorbed by the canopy and the latent heat lost from the canopy. Below the line the canopy is warmer than the air, above the line the canopy is cooler than the air. Compare this with the figure of canopy temperature above. At a relative humidity of 90%, both the blue lines cross the respective reference lines at ca. $110$$110$ W/m². This verifies the consistency of eqs. $\text{27}$$\ref{eq_t_plant}$ and $\text{33}$$\ref{eq_transpiration}$. The radiative energy absorbed by the canopy ends up in latent heat bound in the air and in the heating of the canopy itself.

The convective model of heat transfer does not match the complex physiology and physics of the canopy layer; all heat fluxes are taken care of by eqs. $\text{27}$$\ref{eq_t_plant}$ and $\text{33}$$\ref{eq_transpiration}$. Therefore the canopy layer is given a $U$$U$-value of zero to make it fit into the stack of generic layers.

Heat pipe layer

Heat pipes are installed to heat the greenhouse by convection and long-wave radiation. The transmissivity of the heat pipe layer is ${\tau }_{pipe}={\tau }_{pipe}^{\prime }=1$$\tau_{pipe}=\tau_{pipe}'=1$ for all wave lengths.

The drop in temperature ($\mathrm{\Delta }{T}_{pipe}$$\Delta T_{pipe}$; K) from the inlet to the outlet is modelled by the empirical equation

where $k$$k$ and $b$$b$ are parameters calibrated to the greenhouse, $\mathrm{\Delta }{t}_{pipe}$$\Delta t_{pipe}$ (s) is the transit time of pipe water, $T_{pipe}^{inlet}$$T_{pipe}^{inlet}$ (°C) is the water temperature at the pipe inlet, and $T_{in}$$T_{in}$ (°C) is the greenhouse air temperature. If the inlet temperature is not above the greenhouse temperature then $\mathrm{\Delta }{T}_{pipe}=0$$\Delta T_{pipe}=0$.

The transit time $\mathrm{\Delta }{t}_{pipe}$$\Delta t_{pipe}$ (min) is related to the flow rate ${\dot{v}}_{pipe}$$\dot{v}_{pipe}$ (m³/h) and pipe volume $V_{pipe}$$V_{pipe}$ (m³) as

As an example, a pipe with an inner diameter of 30 mm installed at a density of 2 m/m² in a greenhouse area of 10 000 m² holds the volume,

which at a flow rate of ${\dot{v}}_{pipe}=20$$\dot{v}_{pipe}=20$ m³/h gives a transit time of

The energy lost from the heat pipe $E_{pipe}^{tot}$$E_{pipe}^{tot}$ (W/m² ground) is related to the temperature drop through the pipe $\mathrm{\Delta }{T}_{pipe}$$\Delta T_{pipe}$ and the flow rate ${\dot{v}}_{pipe}$$\dot{v}_{pipe}$,